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Abstract 



A hybrid scheme obtained by combining 3DVar with the Assimilation in the Unsta- 



> 

(N 

O 
00 
O 

> 

X. 

ble Subspace (3DVar-AUS) is tested in a QG model, under perfect model conditions, 
with a fixed observational network, with and without observational noise. The AUS 
scheme, originally formulated to assimilate adaptive observations, is used here to as- 
similate the fixed observations that are found in the region of local maxima of BDAS 
vectors (Bred vectors subject to assimilation), while the remaining observations are 
assimilated by 3DVar. The performance of the hybrid scheme is compared with that of 
3DVar and of an EnKF. The improvement gained by 3DVar-AUS and the EnKF with 
respect to 3DVar alone is similar in the present model and observational configuration, 
while 3DVar-AUS outperforms the EnKF during the forecast stage. The 3DVar-AUS 



algorithm is easy to implement and the results obtained in the idealized conditions of 
this study encourage further investigation toward an implementation in more realistic 
contexts. 

1 Introduction 

Data assimilation in meteorology and oceanography is experiencing a rapid phase of devel- 
opment, with flourishing of theoretical studies and applications. Traditionally the regular, 
globally available measurements network has been employed in the process of combining 
data with model forecasts to obtain the initial condition for model integrations. Recently 
there has been also a growing interest in adaptive observational systems ([22 1) and satellite 
data thinning (|26j). 

While the merits of the relatively new ensemble methods are being compared with those 
of traditional variational methods, the potential benefit of the use of adaptive observations, 
added to the standard network, is being tested ([32]). 

The present is a follow-up study on the testing of an assimilation scheme developed 
by the authors ([3S1 SHI [HH1 E!) and referred to as Assimilation in the Unstable Subspace 
(AUS). The basic paradigm of AUS is to track and control the analysis- forecast cycle 
instabilities. As a consequence, AUS has found its most natural applications in an adaptive 
observation context. Making use of just a few additional observations, properly located at 
each observing time in order to detect the analysis-forecast cycle instabilities estimated by 
Breeding on the Data Assimilation System (BDAS) it was possible to drastically reduce 
the forecast error. 

Previous studies have in fact proven the AUS-BDAS system effectiveness in a hierar- 
chy of dynamical systems characterized by different features and complexity, such as the 
Lorenz 40- variable model ([23]). the atmospheric quasigeostrophic model used here (|_29|) 
and a primitive equation ocean model ([!]). These results support, on the one hand, that 
BDAS, in view of its ability of tracking the instabilities of the assimilation system, provides 



a feasible targeting strategy; on the other hand that AUS is an efficient, computationally 
affordable, relatively easy to implement method to assimilate the adaptively located obser- 
vations. Moreover, they show that, even dealing with high-dimensional system, an efficient 
error control and an accurate state estimate can be obtained by monitoring only a reduced 
number of unstable directions and by properly using a limited number of observations. 

By allowing ad hoc observations of the estimated instabilities, the deployment of an 
adaptive network naturally enhances the efficiency of AUS. In fact, targeting the spatial 
unstable structures improves the estimate of the amplitude of the unstable components of 
the error, thus enhancing the stabilizing effect. One natural question is therefore how AUS- 
BDAS performs with a fixed, standard, observational network. As long as the spatial and 
temporal observational coverage of the standard network is dense enough, the instabilities 
can be detected as they move throughout the domain, although a reduced efficiency may 
be expected with respect to the adaptive observations case. 

|39j investigated the use of fixed "satellite" surface elevation observations in a primitive 
equation oceanic model and pointed out the weakness of such an observational network. 
A stable error reduction could not be achieved by using only surface observations since 
dominating instabilities developed in deeper model layers before they could be detected 
at the surface. Instead, an efficient error reduction and a drastic improvement of the 
assimilation performance were obtained by means of vertical profiles adaptively located 
by BDAS and assimilated by AUS. 

We investigate the best use of the AUS-BDAS method in an atmospheric context where 
only a network of fixed soundings is available. 

Errors, in an operational data assimilation cycle, are of various origin. Even if it is 
assumed that errors in the description of the system dynamics (the model error) are not 
present, noise is introduced at different stages of the observing and assimilation procedure. 
It includes, among others, the measurement noise (coming also, for instance, by deficiencies 
in the retrieval procedures of indirect measurements), representativity errors and errors 



due to imperfections in the assimilation procedure itself. Moreover, if the instabilities 
cannot be controlled by adequate (ad hoc) observations, nonlinearities become important. 
Thus, in the presence of large unstructured errors due to observation and system noise 
and to strong nonlinearities that can be expected to be important in a fixed observational 
setting, a failure of the unstable subspace paradigm of AUS is expected. These arguments 
suggest and motivate the combined use of AUS and a stationary assimilation method. 

In the present study we propose an algorithm composed of a least-square based data 
assimilation scheme, the 3DVar, in combination with AUS, hereafter referred to as 3DVar- 
AUS. The added value of AUS over a 3DVar analysis cycle is investigated with a fixed 
observational network chosen to simulate a realistic observation coverage, and observation 
system simulation experiments are performed with both perfect and noisy observations. 

Given the properties of the available observing network, a practical relevant question 
is how AUS compares with other assimilation methods. For comparison, an ensemble 
Kalman filter (EnKF) ( [121 H~3] ) , has been implemented and optimized for the model and 
observational network used here. The EnKF is considered among the most promising 
ensemble-based sequential data assimilation schemes and therefore is an ideal candidate 
for the comparison with the proposed 3DVar-AUS scheme. 

This article is organized as follows. In Section 2 the details of the model and of 
the observational network are given. In Section 3 the three data assimilation schemes, 
objects of the comparison, are briefly described along with the specific implementation 
choices adopted here; particular emphasis is given to the proposed 3DVar-AUS. Results 
are described in Section 4, while the final summary and discussion can be found in Section 
5. 

2 Model and observation network configuration 

All the experiments described in this study are performed by making use of an atmospheric 
numerical model based on the quasigeostrophic equations in a periodic channel ([29]). The 



model features mid-latitude large-scale flows and the channel is centered at 45°N; its length 
is 16000 km while its width is 8000 km (approximately 180° of longitude by 70° of latitude). 
At the resolution of 250 Km used here, the model fields are discretized at 64 x 33 horizontal 
gridpoints on 5 vertical inner levels where potential vorticity, PV, is defined, plus 2 (top 
and bottom) levels where potential temperature, PT, is defined, so that the model state 
vector has dimension / = 14784. The model solution is forced by relaxation to a zonal 
mean state with constant stratification and damped by a V 4 horizontal diffusion and by an 
Ekman pumping at the surface. The time step for integration is approximately 30 minutes. 
At the resolution used, the model possesses 24 positive Lyapunov exponents, the leading 
one corresponds to a doubling time of 2.2 days. This model provides a qualitatively good 
representation of the true atmospheric mid-latitude dynamics while being simple enough to 
make long runs and statistical analysis feasible; it has been already used in some previous 
works on AUS ([HUE]) as well as in a number of other data assimilation studies ([25ll8|l9]). 
An analysis of its dynamical properties can be found in |30j. 

The forward model integration evolves from the analysis state at a given time t k to the 
forecast state at time tk+i = t^ + t where k = 0, 1, 2..., and r is the assimilation interval 
of a sequential assimilation scheme: 

x/(t fc+1 ) = M( X a (t k )), (1) 

where x-^ and x a indicate the forecast and analysis states respectively. 

The observing network consists of a fixed set of observations simulating synoptic pro- 
files located at model gridpoints and measuring model variables, (potential vorticity, PV, 
at the five inner levels and potential temperature, PT, at top and bottom), at each level. 

In a previous work with the same model ([7]), AUS was used to assimilate a single 
adaptive observation in a large data-void area, and was combined with a 3DVar that 
assimilated all the fixed observations over a data-rich area. 

In the present study, we intentionally work with a fixed network of observations to 
investigate the potential usefulness of AUS in the absence of adaptive observations that 



naturally enhance its efficiency. The observational network used here (Fig. [T]) alternates 
data-rich/data-void areas as it is typical of a real observation network: the distribution of 
M = 125 soundings is obtained through a random selection procedure. 

Observing system simulation experiments are performed: a given trajectory, solution of 
the model equations, is taken to represent the true atmospheric evolution. At each analysis 
time, every r = 12 hours, observations are taken by sampling the "true trajectory" at the 
observation locations. Following the widely used notation (|_20j), observation values are 
stored as components of the observation vector y°: 

y°{t k ) = H(yL(t k )) + e(t k ) (2) 

where e(tfc) represents the additive observation error, assumed to be white in time, Gaus- 
sian distributed with covariance R and TC is the observation operator, estimating the 
observed variables from the model state. Since the observations are located at model 
gridpoints and observe model variables, the observation operator is inherently linear and 
the notation H is used hereafter. 

Observation errors are assumed to be correlated in the vertical direction only and 
uncorrelated for different soundings. The observational Root-Mean-Square (RMS) error 
at each level is set to 10 % of the system's natural variability while, as in Morss (1999), 
the covariance between different levels is expressed in terms of the respective variances 
and of a vertical correlation function (Eq. 3.19 in [I]), depending on the vertical distance 
only. With these assumptions the observation error covariance matrix R takes the simple 
form of a diagonal block matrix, with a 7 x 7 square matrix in each of the M diagonal 
blocks. 

Following [23] and [17], observation errors are randomly generated, and added to the 
true values, by sampling a Gaussian distribution AA(0, R) with zero mean and consistent 
with the assumed observation error statistics. 



3 Data Assimilation Algorithms 

The analysis state at an arbitrary time tj- is obtained as a linear combination of a forecast 
state, taken as a background field, with the observations: 

x a = (I - KH) x f + Ky° = x / + Kd (3) 

where d = y° — Hx' is the innovation vector, while the expression of the gain matrix 
K that minimizes the analysis error variance is the Kalman gain ([21J): 

K = P^H 1 [HP / H T + R.1 ' (4) 

where P^ represents the forecast error covariance matrix. In Eqs. f£5§ and (|3|) all 
vectors and matrices refer to the same time tk, and T indicates matrix transposition. 

3.1 3-Dimensional Variational Assimilation 

The three dimensional variational algorithm, 3DVar, used here has been developed on the 
basis of that described in [23]. The forecast error covariance is estimated by the stationary 
matrix Pl DVar = B (background error covariance), so that the analysis update is (see e.g 
[33]): 

x a = x. f + BH T [HBH T + R] _1 d (5) 

corresponding to the minimum of the 3D-Var objective function for a linear observation 
operator H. The matrix B, that is kept constant throughout the analysis cycle, was 
statistically estimated by [23] and has been multiplied here by a scalar coefficient in order to 
optimize it for the present observational setting. The optimal value, chosen by minimizing 
the analysis error with the present network of 125 noisy profiles (Fig.l) is 0.97. 

Equation © is solved globally by using a conjugate residual solver algorithm ([23]). 
The initial conditions for all the experiments described in the text is the final state of a 
1-year-long 3DVar analysis cycle initialized by randomly perturbing the truth state. 



3.2 Ensemble Kalman Filter 

The ensemble Kalman Filter, EnKF, used here is based on that described in Descamps and 
Talagrand (2007) and used for the study and intercomparison of ensemble initialization 
techniques. It basically follows the approach given by |131I14|. 

Let {x^(tfc)} and {x.f(tj-)}, with i = 1,2, .JV, be the ensemble of N forecast and 
analysis states at times £&, while x^(tfc) and x a (tfc) represent their respective means. The 
forecast and analysis deviations from the mean, 5x.? ,a = x.f' a — 5t-'' a are stored as columns 
of the I x N matrices X-" and X a , respectively. The ensemble-based forecast and analysis 
error covariance matrices are then defined as: 

HnKF = aF37 X/X/T P EnKF = J^J^^ (6) 

At each assimilation time tk, a set of perturbed observations is generated, and then 
used for the ensemble analysis update, by randomly perturbing the actual observation 
vector: 

yt = y° + ei i = l,2,...,N, (7) 

where £{ are independent realizations of the observational noise 7V(0, R). 
At the analysis time the ensemble of forecast states is updated: 

x? = (I - K EnKF H)x{ + K £n ^ F y 4 ° * = 1, 2, ..., N, (8) 

where the gain matrix K.EnKF is obtained by setting P-^ = P EnKF in Eq. (jlj). The 
ensemble of analysis states is then used to initialize N full nonlinear model forecasts. 

In all the experiments described here, N = 30 members are used and the very first 
analysis ensemble is generated by adding random noise to the reference initial condition. 
This noise is drawn from a Gaussian distribution with zero mean and variance equal to 
the average analysis error variance of the 3DVar experiment. 

In order to avoid ensemble collapse, prevent filter divergence and optimize the EnKF 
performance, following [TU], forecast error covariance inflation ([2]) and localization ([ 



have been adopted. The covariance inflation is obtained by multiplying the matrix P^ by 
a scalar coefficient before using it in the analysis: 

PL^a + a)^^, (9) 

where a is the (tunable) inflation factor. As in Houtekamer and Mitchell (2001), the 5th 
order pjy function is used to localize the ^euKF' the ma i n parameter being the zero- 
correlation distance, do- 
Both a and do have been optimized, over a 60 days analysis cycle, for the specific model 
and noisy observational network used here. The best results are obtained with do = 3000 
Km and a 7 % of inflation (a = 0.07). Figure [2] shows the normalized RMS analysis error 
as a function of the inflation factor (fixed zero correlation distance do = 3000 Km), and 
as a function of the zero correlation distance (fixed inflation factor a = 0.07) in the inset. 

3.3 Assimilation in the Unstable Subspace (AUS) combined with 3DVar 

3.3.1 Assimilation in the Unstable Subspace 

The mathematical formulation of AUS is briefly reported in the following, along with the 
description of the setup of its application in the present model and fixed observations 
configuration. A full description of its mathematical and theoretical foundations can be 
found in |38j,|39| and [7\. 

The AUS method is basically aimed at tracking, and controlling, the instabilities of 
the (observationally forced) analysis- forecast cyclic system in view of their prominent role 
in the error evolution and their impact on the overall performance of the data assimilation 
scheme. To this end, when observations are available, an estimate of the analysis-forecast 
cycle unstable subspace is used to confine the analysis correction within this subspace, 
thus maximizing the effect of the assimilation where it is expected to be more necessary. 

By combining Eqs. (1) and ([3]), the evolution equation for the analysis state, the 



analysis-forecast cycle, is obtained: 

x a (t fc+1 ) = (I - KH)A4(x a (t fe )) + Ky°(t fc+1 ), (10) 

Apart from the presence of a linearized observation operator H, Eq. (I10p is the governing 
equation of most sequential data assimilation schemes. 

The analysis-forecast cycle perturbation equations are given by: 

5x a (t k+1 ) = (I - KH)M5x a (t fe ) (11) 

where M, the Jacobian matrix of M, represents the tangent linear model evolution. In 
principle, from Eq. (|11|) . it is possible to estimate the unstable manifold of the forecast- 
analysis system. 

After storing the unstable vectors as the columns of a I x N matrix E, and confining 
the analysis increment in the unstable subspace, the AUS analysis becomes: 

x a = x' + K At/5 d, (12) 

where: 

K AUS = Er(HE) T [R + (HE)r(HEf] _1 , (13) 

or, equivalently: 

K AUS = E [(HEfR-^HE) + T" 1 ]" 1 (HE) T R~ 1 , (14) 

where, as before, all terms refer to the same time tf. and where T is a N x N positive 
definite matrix, representing the forecast error covariance matrix in the unstable subspace: 

Hus = ErET (15) 

Suppose that a single mode e is detected by a set of observations, then the AUS 

analysis reads: 

, (He) T R- 1 d . , 

X = X/ + /TT NTT, 1 /TT X 9 e !6 

(He) T R- 1 (He) +7-2 v ; 

where 7 2 , the forecast error variance in the e direction, can be estimated from innova- 
tions. 

10 



From Eq. (ll6p we see that the analysis increment vector has the direction of e, and 
the amplitude that best-fits the observations. This means that, in physical space, the 
difference between the analysis and the forecast state has the three-dimensional structure 
of the unstable perturbation. 

In practice, the unstable subspace of the data-assimilation system is estimated by 
means of an extension of the Breeding technique ( [35} I36j). known as Breeding on the 
Data Assimilation System (BDAS) (|38j), which naturally incorporates the observational 
forcing in the perturbations dynamics. 

Equation (TTTj) is at the base of BDAS although in practice, as in the standard breeding, 
the full nonlinear model M. is used instead of M to evolve initial random perturbations 
which are then kept small by repeatedly scaling them down to their initial amplitude. At 
each assimilation time, Eq. (|12p is used to update both the control reference trajectory, 
solution of (1), and the set of the BDAS perturbed trajectories. 

This choice is motivated by the interest in testing and developing techniques feasible 
for realistic contexts where coding, implementing and keeping updated a tangent linear 
model may require some effort (for details on the implementation of BDAS see |38^ HO] 
and Sect. 3.3.2). 

Furthermore, it is impractical to estimate the whole unstable subspace; this would 
require either a recursive orthonormalization to be applied to the set of BDAS vectors, 
according to the standard technique of [3] or the use of computationally demanding tech- 
niques such as those described [37] and by [35] , which allow the estimation of the Lyapunov 
vectors. Instead, in this as in previous applications of AUS, a refresh procedure is intro- 
duced in the breeding cycle in order to systematically explore the unstable subspace of 
the system. In fact, as discussed in [38J and [39], the more accurately the forecast error 
projection on a particular unstable direction is estimated, the more dominant will become 
the error projection on the complementary subspace. 

Once a breeding cycle has been established, the BDAS modes and the forecast error 
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will approximately have a structure given by a different linear combinations of independent 
unstable vectors of the underlying dynamics. 

Local structures appearing in the global BDAS modes will be positively or negatively 
correlated with similar structures in the forecast error. It is therefore necessary to localize 
these structures before entering them in an assimilation expression of type (16). In addi- 
tion, the localized structure needs to be observed, i.e. He should be non-zero and large 
enough to provide a reliable estimate of the analysis increment, whose sign is correctly de- 
termined by the innovation^. It is important to notice that, by applying only a horizontal 
localizing function, the vertical structure of the perturbation is preserved in the definition 
of the analysis increment. 

While EnKF-like methods are based on an estimate of the forecast error covariance 
built from an ensemble of trajectories, the identification of unstable directions of the given 
control solution, which provide the three-dimensional structure of the analysis increment 
whose sign is determined by the innovation, is the distinctive feature of AUS. 

EnKFs are based on the traditional Kalman filter equations and use a Monte Carlo 
approach to estimate the forecast error covariance matrix entering the analysis update, 
therefore they have to face the consequences of sampling errors with or without perturbed 
observations ([43]) and filter divergence. As a consequence, several ensemble filters for- 
mulations have been proposed to overcome these difficulties (for a review see PQ, [13] . 
|34J). Similarly, the AUS approach, that relies on the estimate of the forecast error pro- 
jection on the unstable subspace, faces the practical difficulties of accurately estimating 
the unstable directions and the amplitude of the forecast error on each of them. While the 
theoretical premises of the two approaches are different but clear, when it comes to the 
implementation all schemes are subject to several practical choices that make quantitative 

comparisons rather subjective. 

1 By using R _1 as metrics, the norm of He is (He) R" 1 (He): the components of e detected by 
different observations are weighted by the observational accuracy. Noisy observations can then be effective 
if the e component they detect is large enough. 
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An adaptive observations approach was followed by [39] and by [7] in an oceanic and 
atmospheric model applications respectively. In fact, the AUS approach of tracking and 
controlling the data assimilation system instabilities is most efficient when observations 
are indeed available at the moment and in the region where an instability develops. In 
other words, if the set of unstable vectors E can effectively be detected (that is to say if 
the components of HE are large enough), then E can be used in Eq. (1121) to update the 
background state. When, as in the present application, only fixed synoptic observations 
are available, the efficiency of the AUS scheme is expected to be determined by the spacing 
and frequency of the observational network in relation to the typical spatial patterns and 
growing times of the system's instabilities. 

3.3.2 Implementation of 3DVar-AUS 

Several approaches have been proposed to combine 3DVar with ensemble-based filters. 
These hybrid Ensemble/3DVar analysis schemes introduce information on flow-dependent 
instabilities by properly weighting the static covariance with the ensemble estimated fore- 
cast error covariance (|16]l41j). [llj investigated the robustness of hybrid schemes to model 
error: the model error modifies the optimal values of the weights to be given to the ensem- 
ble based and static covariance matrices, significantly reducing the weights corresponding 
to the former. 

In the present study, we introduce a hybrid scheme that uses the flow-dependent co- 
variance and the 3DVar covariance separately to assimilate different observations. The 
AUS assimilation is applied first to reduce the forecast error component in the model's 
unstable subspace estimated by the BDAS modes, using those observations that are able 
to detect them. The residual forecast error is assumed to be unstructured and a static 
3DVar covariance is used to assimilate the remaining observations. Possibly, error as- 
sociated with unstable structures that were not present in the BDAS modes, or whose 
amplitude could not be reliably estimated with the available observations may be still 
present in the analysis and subsequent forecast: hopefully they will be accounted for at 
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successive assimilation times. 

The 3DVar-AUS forecast-assimilation and breeding cycle is implemented through the 
following recursive steps: 

1. add N (= 1) random perturbations to the control analysis state (store the new set 
of perturbed states in place of those discarded; see step 9); 

2. perform forecast integration of the control analysis and all perturbed states. The 
total number (= 12) of perturbed trajectories evolved is equal to the breeding time 
(6 days) divided by the assimilation interval (12 hrs) times N (= 1); 

3. determine the N(= 1) global BDAS modes to be used in the following step, as the 
difference between the perturbed trajectories that have completed a forecast- analysis 
breeding cycle (6 days) and the control forecast; 

4. select from the N(= 1) global BDAS vectors the (N°, see below) local structures 
that are detected by observations; 

5. estimate the forecast error variance, 7 2 , from recent innovations; 

6. apply AUS analysis, using the selected structures with relative observations, and 
making use of 7 2 (step 5), to the control and to the full set of (= 12) perturbed 
trajectories; 

7. store the innovations to be used in step 5; 

8. apply 3DVar, using the remaining observations, to the control and to the full set of 
(= 12) perturbed trajectories; 

9. discard the trajectories relative to the N(= 1) BDAS modes used in steps 4-5, and 
rescale the remaining perturbations; 

10. go back to step 1. 
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The values given in parenthesis refer to the specific values adopted in the present 
application. 

During the first 6 days of experiment, the control and all the perturbed trajectories 
are subject to only the 3DVar analysis update to assimilate all the observations; the AUS 
analysis cycle begins at day 7. 

Thus, in all the experiments described below, a single N = 1 BDAS mode is used at 
each assimilation time while the breeding time interval is set to 6 days. This latter value 
was chosen by tuning to minimize the average analysis error (see also [7J). The breeding 
time is related to number and growth rate of the unstable directions of the reference 
trajectory, solution of the assimilation cycle. In practice the breeding time has to be long 
enough for the perturbations to acquire the structure of the local instabilities of the flow. 

At the analysis times, instead of directly using the global BDAS modes e in Eq. (16), 
N local isolated structures are extracted from e. Once the unstable structures are selected, 
we identify the observations, if any, able to detect the N local structures; this defines a set 
of N° < N "observable" structures which are finally used in the analysis update. Details 
on steps 4-6 are given in the following: 

• (a) Selection and localizations of structures. The localization is made by repeat- 
ing the point-by-point multiplication of an horizontal Gaussian function, /(x) = 
g||x— xmoaell /d ^ being the characteristic distance), centered on the local maxima, 
Xmax, of each BDAS mode e. The process starts from the absolute maxima; sub- 
sequent maxima are searched at a distance larger than twice d. The N extracted 
structures are therefore assumed to be uncorrelated. We point out that the local- 
ization is made in the horizontal direction only, so that the vertical structure of the 
local maxima is preserved. 

• (b) Identify available observations and observable structures. The observations se- 
lection is made in two steps: all the structures with no observations within an 
horizontal (1,1) gridpoints box, centered on their maximum, are discarded; the am- 
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plitude of the structure at the observation locations must be larger than a fixed 
ratio /3 of its maximum amplitude. At the end of these steps, N° < N "observable" 
structures, e°, are obtained, each one detected by Mi observations at the same level, 
i = 1,...,N°. For each of the N° observable structures, define one innovation vec- 
tor, dj, one observation error covariance, the Mi x M{ matrix Rj, one observation 
operator, the M{ x / matrix Hj and one scalar 7?. 

(c) AUS analysis update. If N° 7^ 0, that is at least one observable structure e° is 
identified, the AUS analysis update is performed, otherwise all the observations are 
assimilated with 3DVar. The analysis is made according to Eq. (16), using the N° 
structures e° sequentially: 

x a = xf + V (H i e°) T (R)- 1 d t 

+ ^(H l e^(R l )-i(H l e°) + ( 7 ,)- 2e4 {U) 

Expression (fT7|) is the AUS analysis equation in the case of N° isolated structures 
each one observed by Mi observations whose associated error covariance matrices 
are Rj. 

Given the characteristics of the observational network used here and in view of 
the selection procedure (point (b)) , each set of Mi observations associated to the 
observable structure e° consists of scalar measurements at the same level (where the 
structure maximum is located). As a consequence, the observations are uncorrelated 
and all have the same variance a 2 : R» = cr 2 Ij, where Ij is the identity matrix of order 
Mi. By making use of these hypotheses, after rearranging, Eq. (|17|) becomes: 



U a2+ ^(H^HH^) (H^Hh^) < l 5j 

Equation (|18|) is the AUS analysis update used in the experiments described here- 
after. Each of the scalars 7 4 2 , representing the variance of the forecast error along 
the structure e°, needs to be estimated. 
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In practice the terms 7? (Hje°) (Hje°) are estimated statistically using the inno- 
vations, according to the approach introduced in [7j: details can be found in the 
Appendix. 

Once the AUS analysis is completed, the analysis field is used as the background in the 
3DVar analysis to assimilate the remaining M — X^j=i ^i observations (step 8). 

As mentioned above, a refresh procedure (steps 1 and 9) is used in the implementation 
of BDAS ([38]). At each assimilation time, after the analysis update, the set of BDAS 
modes is discarded after being used in the assimilation. A new one is introduced, which 
starts to undergo a new breeding procedure: after completion of its breeding period it will 
be used in an assimilation step. Although this procedure increases the computational cost 
of BDAS (with respect to evolving the same set of perturbed states), it is very beneficial to 
efficiently span the trajectory's unstable subspace and to provide a reliable estimate of the 
data assimilation system instabilities. The refresh was also used in [7] in the context of the 
same QG model used here, while a trade-off procedure between keeping and discarding the 
iV(=6 in that case) BDAS modes altogether was used in [39] in the context of a primitive 
equation ocean model. 

In the experiments described here, the characteristic distance is chosen, after optimiza- 
tion, to be d = 1500 Km. Also, a limit to the total number of the extracted assimilating 
structures is set, N = 20. Furthermore, the size of the searching box I and the coefficient 
[3 have been optimized by tuning, and are equal to I = 7 gridpoints (equivalent to an area 
of 1500 Km 2 ) and [3 = 0.6 

A practical aspect which has to be taken into account when comparing 3DVar-AUS and 
EnKF is the computational cost. Both schemes, the hybrid 3DVar-AUS and the EnKF, 
require high computational power determined mainly by the number of evolving trajec- 
tories (perturbations of control for the 3DVar-AUS; ensemble members for the EnKF). 
At the analysis step, the computational most demanding part in the EnKF are the large 
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matrix multiplication and inversion needed to compute the gain matrix; in 3DVar-AUS, 
the computational cost is mainly due to the 3DVar analysis update which has to be re- 
peated for all the simultaneously evolved perturbed trajectories. A detailed comparison 
of the computational cost of the two algorithms needs to be done case by case, since the 
implementation choices are expected to depend on the system under consideration and on 
the available observation network. The computational cost, in terms of computing time, 
in the present application are given at the end of the next section. 

4 Results 

4.1 Experiments with perfect observations 

We first compare the performance of the algorithms in a perfect observation setting. 

Clearly real observations are not perfect. Anyhow, the idealized perfect observations 
setting is related to the problem of observability of the assimilation system and provides 
insight on the ability of various methods to track the instabilities of the system. The use of 
perfect observations with AUS allows for an accurate evaluation of the analysis increment 
and potentially reduces to zero the analysis error projection along the unstable structures 
used in the analysis update (18); the relation between the observability condition and 
AUS was discussed by |38j. Furthermore, the perturbed observation case provides an 
upper limit of performance of the assimilation schemes under comparison. 

Perfect observations are obtained by applying the observation operator to the true 
state; e = in Eq. ([2]) and £j = in Eq. ([7|). The observation error covariance matrices 
are set to zero in Eqs. © and ([H]) and a 2 = in (|18|) . Under these conditions all EnKF 
members assimilate the same (unperturbed) observations y°; the parameters a and efo °f 
the EnKF experiments are set to the same values optimized in the noisy case, a = 0.07 
and d = 3000 Km. 

Figure [3] shows the normalized RMS analysis error, as a function of time, for 3DVar, 
3DVar-AUS and EnKF over 2 years. The 3DVar analysis error undergoes fluctuations 
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with large spikes, presumably related to an unreliable representation of the actual forecast 
error in the stationary estimate B of the forecast error covariance matrix. 

The average RMS analysis error, excluding the first 100 days, is about 6 % of the 
system's natural variability. The behaviour of 3DVar-AUS and EnKF is impressive: the 
average RMS analysis error, excluding the first 100 days, is about 1.6 x 10~ 6 (0.00016 %) 
and 5.0 x 10 -7 (0.00005 %) of the system's natural variability, respectively. At the very 
beginning of the assimilation cycle, the EnKF solution drifts away from the truth, probably 
because the P EnKF has not yet acquired the dynamic consistency, then, within about 50 
days, the analysis error rapidly decreases and remains confined to very low values for all 
the subsequent period. Analogously the 3DVar-AUS experiment takes a transient time 
before error reduction; after about 80 days the 3DVar-AUS analysis error drops to very 
low values until the end of the experiment. 

Figure S] shows the current BDAS mode superimposed to the actual forecast error, at 6 
instants during the 3DVar-AUS perfect observations experiment. Although the coincidence 
is not systematic, the BDAS mode ability to capture the relevant components of the 
forecast error appears evident. The actual forecast error patterns appear composed of 
spatially distributed structures, often in a dipole or quadrupole configuration, separated by 
wide, relatively homogeneous, low error regions. Such pronounced and localized structures, 
presumably related to the instabilities of the underlying dynamics, are also frequently 
found in the corresponding BDAS mode. These are the instabilities we want to track and 
control. If they are not eliminated by the assimilation at one analysis time, they may still 
migrate to observed areas and, possibly, be detected at a later time. From Fig. [5] we also 
see that the spatial scale of these structures does not exceed 3000 Km (11 gridpoints): 
this supports the use of the Gaussian masking function with d = 1500 Km. 

The local character of the atmospheric instabilities are described in [28] and have 
been efficiently exploited in AUS as well as in most of the ensemble prediction and data 
assimilation applications, as for instance the Local Ensemble Kalman Filter (|27 |. I3"T | 19]) 
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and the Local Ensemble Transform Kalman Filter ( [19] ). 

4.2 Experiments with noisy observations 

Although very instructive and useful in providing insight into possible advantages and 
drawbacks of a given algorithm, the use of perfect observation is not realistic: in real 
applications, the assimilation of noisy observations introduces unstructured random errors. 
Therefore, we now turn to noisy observations experiments; observations are generated 
according to the procedure described in Sect.EJ A set of observation realizations, Eq. ((7|), 
is used to update the ensemble forecast in the EnKF by means of Eq. (jHJ) , while Eq. (fT8l) 
is used for the AUS assimilation in the 3DVar-AUS experiment. 

Figure [5] shows the normalized RMS analysis error, as a function of time, relative 
to the 2 years trajectory obtained by the assimilation of noisy observations with 3DVar, 
3DVar-AUS and EnKF. The 3DVar RMS error, excluding the first 100 days, is now about 
22 % of the system's natural variability, more than twice the RMS observation error, and 
large spikes are still present. Remarkably, both 3DVar-AUS and EnKF have an average 
RMS analysis error of about 7 % of the natural system's variability, below the RMS 
observational error (10 %, Sect. [2]). 

Reducing the number of ensemble members to 25 and 20 (optimized values do = 3000 
Km and a = 0.09 in both cases), the EnKF RMS analysis error is equal to 10.08 % and 
14.21 % respectively. 

Notice that besides having comparable time mean values, the 3DVar-AUS and EnKF 
RMS analysis error shows highly correlated fluctuations which can be interpreted as in- 
duced by the same flow-dependent instabilities of the underlying dynamics. The large 
improvement of the 3D Var-AUS over the standard 3DVar appears impressive considering 
that a single BDAS mode and very few observations are used by AUS at each assimilation 
time. 

During the 2 years of simulated time, the AUS analysis was performed 89 % of the 
total assimilation times: the analysis update was based, on average, on N° = 9 observable 
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structures e° each one detected by (most of the times) Mj = 1 observation. 

These findings confirm the authors' claim that a limited number of observations and 
unstable directions, if properly used, allow to improve upon the 3DVar performance and 
make the hybrid scheme performance comparable to that of the EnKF. 

Figures 6 and 8 show the 3DVar (left column) and the EnKF (right column) analysis 
increment superimposed to the 12 hrs actual forecast error, in a sequence of three assimi- 
lation times starting at day 204, 00UTC. For the same times, Figs. 7 and 9 show, for the 
3DVar-AUS experiment, the increment obtained by AUS only (left column) and 3DVar 
(right column) superimposed to the 12 hours forecast error and to the background error 
after the AUS update. 

In all panels, black dots indicate the locations of the observations used by the assimi- 
lation algorithm (at the considered assimilation time). Figures 6 and 8 show the top PT 
while the mid level PV is shown in Figs. 7 and 9. 

The 3DVar increments have the typical almost isotropic features coming from the 
simplified assumption in the definition of the stationary matrix B. The EnKF increments, 
on the other hand, reproduce very well the features present in the forecast error: by 
making use of all the available observations and exploiting the accurate description of 
the actual forecast error realized by the 30 members ensemble-based P EnKF-> the EnKF 
provides reliable analysis updates. 

The 3DVar-AUS increments reveal a number of remarkable features. From Fig. 7, we 
see that the number of observations actually used by AUS is very small; 15 scalar obser- 
vations only are used at day 204, 00UTC and 9 observations only at days 204, 12UTC 
and 205, 00UTC. For most of the prominent local maxima in forecast error the AUS anal- 
ysis increment has the proper spatial structure that leads to their reduction. At day 204, 
00UTC, for instance, large error spikes are present along the mid latitude of the channel. 
Two areas are clearly identifiable: one, north-south oriented, between the 5th and the 10th 
longitudinal gridpoints and another, located between the 35th and the 55th longitudinal 
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gridpoints. Both areas are characterized by positive and negative error elongated struc- 
tures. The AUS analysis correction accurately reproduces most of these structures; just 
a small number of observations is sufficient to obtain a reliable and dynamically consis- 
tent analysis. For instance, the elliptic-shaped, westernmost maximum, located between 
longitudinal gridpoints 5 and 10, is accurately accounted for by the localized structure 
extracted from the BDAS mode. However the figure indicates that although the analysis 
increment has the proper shape, residual background error is present as shown in the 
corresponding right panel. Analogously, the strong error maximum at the east side of the 
channel, between longitudinal gridpoints 45 and 55, is reproduced by the AUS increment 
and partially reduced by using a single observation only. The figure also shows a number 
of other forecast error local maxima, some of which captured by AUS. This is the case, for 
instance, of two small dipole-type corrections located between 10 - 15 longitude, 18 - 21 
latitude and between 40 - 45 longitude, 3-9 latitude, each obtained by assimilating two 
observations. 

The need to use the 3DVar after the AUS analysis is highlighted by the right panels 
of the figure. In fact by exploiting the remaining observations 3DVar further reduces the 
error throughout the domain, and not only in the localized areas of large error development 
where AUS is most effective. 

After 12 hours, at day 204, 12UTC, there is no evidence of the large error structure 
which was present in the west side of the channel at the previous analysis time since the 
AUS assimilation efficiently reduced the forecast error in that area, inhibiting further error 
growth. 

The error maximum located between 12 - 20 longitude, 6-16 latitude, "ignored" by 
the assimilation, is still present in the forecast error field. 

The presence of the strong positive error maximum located, as at the previous analysis 
time, in the eastern side of the channel, at mid latitude, suggests that the analysis cor- 
rection at day 204, 00UTC, was not sufficient to eliminate it. Now, at day 204, 12UTC, 
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although its shape is again well captured by the current BDAS mode, the AUS analysis is 
able only to reduce it. Notice that, a very similar pattern (with opposite sign), manifesta- 
tion of the same instability, is present in the EnKF forecast error. Similar considerations 
apply to the last instant shown, day 205, 00UTC. 

For all three algorithms, the mid-level PV increments shown in Figs. 8 and 9 present 
features similar to the PT analysis increments of Figs. 6 and 7. It is most remarkable, 
however, that, except for the single observation at day 204, 12UTC, no observations are 
used by AUS to provide the analysis update at this level. Still, the analysis increments 
strongly correlate with the forecast errors, but the analysis is based on observations located 
at the top level: it is only due to the vertical correlation between the BDAS mode and 
the forecast error that the analysis update is accurate also at levels different from the 
observed ones. The single observation at day 204, 12UTC detects the error maximum on 
the northern side of the channel, between longitudinal gridpoints 21 and 26. 

Figures 7 and 9 illustrate the basic mechanism and the practical conditions necessary 
for the success of AUS: the ability of the observational network, in terms of its spatial 
and temporal distribution, to efficiently detect the unstable structures that grow along 
the trajectory of the assimilation system. They also show the successful use of just a 
single BDAS mode and a small number of observations at each assimilation time. The 
3DVar analysis update, in this hybrid scheme, allows for the use of all the remaining 
observations, distributed throughout the domain to correct errors in regions where AUS 
was not effective. 

The relative contribution of AUS and 3DVAR in the error reduction in the hybrid 
scheme, is illustrated by Fig. 10. The RMS forecast error, AUS analysis error (i.e. the 
background error for 3DVar) and final 3DVar-AUS analysis error, are shown as a function 
of time, for a 50 days period starting at day 175. The plot shows that AUS and 3DVar 
give a comparable contribution to error reduction. This gives further evidence that when 
instabilities are present the few observations used by AUS (approximately 1% of the total 
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available observations) are almost as effective as the remaining ones. In any case, while 
the 3DVar always reduces the corresponding background error of approximately the same 
amount, the effect of AUS is more intermittent. 

Finally, Fig. 11 shows the time and domain averaged analysis and forecast errors as 
a function of the forecast lead time. These are all deterministic forecasts; in particular 
the EnKF forecast is initialized from the EnKF analysis, mean among analysis ensemble 
members. It is interesting to note that although the EnKF and 3DVar-AUS average 
analysis error is comparable, the 12 hours forecast error growth is much more rapid in the 
former case, while it has approximately the same rate afterward. 

With the specific setup choices adopted here, the computing time required to complete 
the 2-year 3DVar-AUS assimilation experiment described in the text, was about 60 % of 
the time required by the EnKF. 

5 Summary and Discussion 

5.1 Summary 

A combined data assimilation scheme, the 3DVar-AUS, based on 3DVar and on the dy- 
namically based AUS algorithm, was presented and discussed here. The test ground of 
the study was the implementation of observation system simulation experiments with a 
synoptic-like network of observations in the context of an atmospheric quasigeostrophic 
model. The performance of this proposed scheme was compared with that of the advanced 
ensemble Kalman filter algorithm. 

According to the formulation given here, 3DVar-AUS is a two-step scheme in which the 
analysis obtained with AUS is used as the background for the 3DVar analysis update. At 
each analysis time, an automatic procedure searched for observable structures, extracted 
from the current BDAS mode, on which the analysis update is based. 

In all the experiments described in the text, just a single global BDAS mode and only 
a very small number of observations were used by AUS at each assimilation time. 
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The system at hand has 24 positive Lyapunov exponents, corresponding to 24 indepen- 
dent diverging directions and inspection of the full spectrum shows that these directions 
have competitive growth rates [6]. These competitive instabilities are present simultane- 
ously in the single global BDAS mode used at a particular analysis time under the form 
of localized structures dominated by local instabilities; these are the structures exploited 
in the AUS assimilation. 

A key ingredient for the success of AUS is clearly the accuracy of the BDAS modes to 
explain relevant components of the actual forecast error. The BDAS procedure is built to 
estimate the instabilities growing along the complete analysis-forecast cycle and implicitly 
embeds all the information about the observational forcing. As indicated by Fig. 0] and, 
indirectly through the analysis increment in Figs. 6 and 7, the correlation between the 
single BDAS mode and the actual forecast error is evident. In particular, the BDAS modes 
reliably reproduce most of the small scales features of the error spatial patterns. 

In the model and observational context used here, the EnKF and the 3DVar-AUS 
showed a comparable performance in terms of accuracy of the state estimate, either with 
perfect or noisy observations. In the perfect observations case, both algorithms are able 
to reduce errors to a very low level demonstrating, on the one hand, an efficient tracking 
of the instabilities (AUS), and, on the other hand, a reliable time-dependent forecast error 
covariance description (EnKF). With noisy observations, the analysis error of 3DVar-AUS 
and of EnKF rapidly decreases, then remains below the observational noise level: the 
improvement over the regular 3DVar alone is dramatic in both cases. 

It is remarkable that, despite the fact that 3DVar-AUS and the EnKF attain a compa- 
rable average analysis error (6.8 % and 6.6 % of the system's natural variability, respec- 
tively) the 12 hours forecast error in the EnKF is larger than in the 3DVar-AUS (9.9 % 
and 7.3 %, respectively) and its 24 hours forecast skill is somewhere between the 36 hours 
and the 48 hours forecast skill of the 3DVar-AUS.. 

This result confirms the efficiency of AUS in reducing the flow instabilities responsible 
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for the error growth even in a fixed observations setting. 

5.2 Conclusions 

The results of the present study lead to the following interpretation: the forecast error 
shows very intense and localized structures that are similar in the three assimilation ex- 
periments and in particular between the 3DVar-AUS and the EnKF: accordingly, a high 
correlation in the time evolution of the RMS error of the two schemes is observed. As 
indicated by these similarities, the local maxima of the forecast error are presumably 
concentrated in regions where dynamical instabilities, compatible with the observational 
forcing, develop along the analysis-forecast cycle. 

With a fixed observational network, observations are assimilated by AUS only when 
they are found in correspondence to the local maxima of the unstable structures. Therefore 
it is necessary to use 3DVar to assimilate the remaining observations. Three dimensional 
error structures, provide the corresponding AUS analysis correction, whose amplitude and 
sign are determined from the innovation: hence its ability to reduce the forecast error also 
at levels other than the observed ones. 

In real world applications, the employment of AUS is expected to give an improvement 
over a 3DVar (or Optimal Interpolation) every time an unstable structure is observed by 
either fixed or adaptive observations. When the latter are available and there is no need 
to wait for the instability to travel into observed regions, the expected improvement in 
using AUS is enhanced. 

Also the EnKF takes advantage of the observations that happen to be (or are inten- 
tionally located) in regions where instabilities develop and the present fixed observations 
experiments prove that the observations are exploited by the two methods with similar 
efficiency. The EnKF obtains the same goal of capturing and controlling the instabilities, 
the different approach of the two methods being exemplified by the following idealized 
example where a single unstable structure is present in the forecast error. In this simple 
case, the EnKF would identify the unstable structure as difference of the members (two 
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members would be sufficient for this purpose) from the ensemble mean. To estimate the 
amplitude of the correction, however, a sufficiently large number of representative mem- 
bers would still be necessary to provide a reliable estimate of the associated error variance. 
In contrast, AUS estimates the direction as difference of a single unstable perturbation 
from the control and estimates the amplitude of the correction from innovation. 

The present result, where EnKF and 3DVar-AUS have a similar analysis performance, 
seems to indicate that both methods were able to efficiently exploit all the available obser- 
vations useful to control the flow-dependent instabilities. A point in favor of the 3DVar- 
AUS scheme is its better performance in the forecast stage. 

Other hybrid schemes have been proposed which combine 3DVar with the ensemble 
based filters ( pH El SD H2] ) • 

Apart from the differences among the different approaches proposed in previous studies 
(see [42] and references therein) two results deserve to be mentioned here. Usually, in these 
hybrid schemes, the forecast error covariance matrix is obtained by combining a static 
covariance with an ensemble based one, through a tunable scalar weight. Improvements 
over a 3DVar conveyed by the flow-dependent covariance are less important when a dense 
observational network is available ( [16] ) . Moreover the presence of model error significantly 
reduces the optimal value of the weight to be given to the flow-dependent covariance (|llj). 

Based on these works and on the results of the present study we draw the following 
conclusion. 

To hybridize an ensemble based scheme or AUS with a static covariance might turn 
to be more convenient when instabilities loose some of their significance, either because 
the observational network is particularly dense ([B]) or if the model error is such that the 
model does not provide an adequate representation of the true unstable modes ( [H] ) . In 
view of the latter consideration, while the hybrid approach proposed in the present study 
may turn out to be useful also in the presence of model error, further work is needed to 
address this problem. 
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Because the performance of different schemes depends upon many different factors, 
including model error and the number and distribution of observations, it remains to 
be seen which method is more efficient in a particular operational setting. Therefore a 
quantitative comparison relevant for operational purposes is out of the scope of the present 
study. 

However, given the ubiquity of the role of instabilities in degrading the analysis accu- 
racy, we believe that AUS may turn out to be useful in those circumstances in which the 
observational network, and in particular an adaptive one, allows their detection. 
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Appendix - Estimate of the forecast error variance along the 
BDAS mode 

In the AUS analysis update Eq. (|18p . each of the scalars 7?, representing the variance of 
the forecast error along the structure e°, need to be estimated. Following the approach 
used in [7] this is done statistically, by making use of the innovations. 
In practice the terms 7?(Hje°) T (Hje°) are estimated as follows: 

i^df^T-* 2 , ^(df,d,) T >a 2 



7/(H i e°) J (H 4 e°) « 



(19) 



i(df,d t ) T , ±(df,d t ) T <a 2 

where (, )t represents the average over an appropriate time interval T and D is a scalar 

coefficient to be tuned. After optimization, these values are set to D = 1.35 and T = 8 

days. 
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Based on Eq. (18), if He approaches zero the analysis correction approaches infinity. 
However because only observations in proximity of maxima in the BDAS mode are chosen 
and in view of the selection procedure described in Sect. 3.3 (point (b)) this circumstance 
is never encountered. 
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Figure 1: The observational network used in the experiments. Each dot represents the 
horizontal location of a vertical sounding profile. 
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Figure 2: Normalized time and domain RMS analysis error for the EnKF experiment as 
a function of the inflation factor (a) by setting do = 3000 Km, and, in the inset, as a 
function of the zero correlation distance do, by setting a = 0.07. The experiments last 
60 days and the average refers to the last 50 days. Errors are normalized with natural 
variability and expressed with potential enstrophy norm. 
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Figure 3: Perfect observations: normalized RMS analysis error as a function of time for 
3DVar (black), 3DVar-AUS (green) and the EnKF (red). Errors are normalized with 
natural variability and expressed with potential enstrophy norm. 
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3DVar-AUS, day 100 00UTC, PV lev 3 



3DVar-AUS, day 200 OOUTC, PV lev 3 
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Figure 4: Perfect observations: mid-level PV background error (shaded) and PV bred 
mode (contour) at days: 100, 200, 300, 400, 500, 600 along the perfect observation 3DVar- 
AUS experiment. 
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Figure 5: Noisy observations: normalized RMS analysis error as a function of time for 
3DVar (black), 3DVar-AUS (green) and the EnKF (red). Errors are normalized with 
natural variability and expressed with potential enstrophy norm. 
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3DVar, day 204 00UTC, PT top 



EnKF (N=30), day 204 OOUTC, PT top 
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EnKF (N=30), day 204 12UTC, PT top 
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Figure 6: Noisy observations experiments: top level potential temperature (PT) 12hr 
forecast error (shaded) and analysis increment (contour) in a sequence of three assimilation 
times starting at day 204, for 3DVar (left cg^imn) and EnKF (right column). Black dots 
indicate locations of the observations used. 
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Figure 7: Noisy observations experiments: top level potential temperature (PT) 12hr 
forecast error (shaded) and analysis increment (contour) in a sequence of three assimila- 
tion times starting at day 204, for 3DVar-AUS. The rights panels show the AUS analysis 
increment superimposed to the 12 hours forecast error; the left panels show the 3DVar 
analysis increment superimposed to the AUS analysis error (which is used as the back- 
ground in the 3DVar analysis update). Black dots indicate locations of the observations 
used. 
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3DVar, day 204 00UTC, PV lev 3 
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Figure 8: The same as Fig. 6 but the mid level potential vorticity (PV) field is shown. 
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3DVar-AUS, day 204 00UTC, PV Iev3-Anlncr AUS 3DVar-AUS, day 204 OOUTC, PV Iev3-Anlncr 3DVar 
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Figure 9: The same as Fig. 7 but the mid level potential vorticity (PV) field is shown. 
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Figure 10: Normalized RMS analysis error as a function of time for the 3DVar-AUS 
experiments. Noisy observations. 
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Figure 11: Noisy observations experiment. Time and domain normalized RMS analysis 
and forecast error as a function of the forecast range: 3DVar (triangles), 3DVar-AUS 
(circles), EnKF (squares). 
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